Center for Turbulence Research 
Proceedings of the Summer Program 1996 

A-priori testing of sub-grid models for chemically 
reacting nonpremixed turbulent shear flows 

By J. Jimenez 1 , A. Lilian 2 , M. M. Rogers 3 AND F. J. Higuera 2 

The /3-assumed-pdf approximation of (Cook & Riley 1994) is tested as a subgrid 
model for the LES computation of nonpremixed turbulent reacting flows, in the 
limit of cold infinitely fast chemistry, for two plane turbulent mixing layers with 
different degrees of intermittency. Excellent results axe obtained for the computa- 
tion of integrals properties such as product mass fraction, and the model is applied 
to other quantities such as powers of the temperature and the pdf of the scalar itself. 
Even in these cases the errors are small enough to be useful in practical applica- 
tions. The analysis is extended to slightly out of equilibrium problems such as the 
generation of radicals, and formulated in terms of the pdf of the scalar gradients. 
It is shown that the conditional gradient distribution is universal in a wide range of 
cases whose limits are established. Within those limits, engineering approximations 
to the radical concentration axe also possible. It is argued that the experiments in 
this paper are essentially in the limit of infinite Reynolds number. 

1. Introduction 

The computation of turbulent reacting flows is an open challenge even after having 
been the focus of intensive work for several decades. The subject of the present note, 
nonpremixed flames with fast chemistry, was one of the first to be tackled, and it 
is somewhat simpler than others, but it still represents a large number of cases of 
theoretical and practical importance. Recent reviews can be found in (Bilger 1989, 
Libby & Williams 1994). 

Our analysis is subject to several simplifications. The diffusion coefficients of 
all the species and of heat are assumed to be identical, k* = k and, although 
not explicitly needed for most of the theoretical arguments, all of our numerical 
experiments are done at Schmidt number Sc = 0.7. Our flows axe incompressible, 
and we assume that any heat released by the reaction is weak enough for the fluid 
density to be unchanged, p = 1. The role of the chemistry is thus passive with 
respect to the flow, although it is modified by it. 

In most cases we assume an irreversible binary reaction 

uaA + vbB — s ► vpP, ( 1 ) 

in a sheax flow between two streams, each of which initially contains either pure A 
or B reactant. 
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Under those circumstances the mass fractions F, of the different species can be 
linearly combined to form a set of conserved scalars which are transported by the 
flow with the common diffusion coefficient k. If in addition the Damkohler number, 
which measures the ratio of the characteristic diffusion and chemical times, is large 
enough, the reaction occurs in a thin flame that can be treated as a surface, and 
the problem reduces to the mixing of a single conserved scalar 


Yao ~Ya + rY B 
Yao + rYg o ’ 


(2) 


called the mixture fraction, which takes values £ = 0 and £ = 1 at the free streams 
(see Williams 1985). Here Y , o is the mass fraction of the i-th species at the appro- 
priate free stream, and 

r = v a WaIvbWb, (3) 

where u, and W, are stoichiometric coefficients and molecular weights. The flame 
is located at the stoichiometric mixture fraction 


“ mw <*> 

and most quantities of interest can be computed as algebraic functions of £, which 
are continuous but have discontinuous derivatives Thus the mass fraction of the 
product Yp is proportional to the triangular function 


/(£) = £/£» if £ < 6, (1 — £)/( 1- W otherwise. (5) 

In modeling turbulent flows we can usually estimate averaged or locally filtered 
values of £, and we would like to have similarly filtered values of mass fractions or 
other quantities, but we are prevented from doing so by the nonlinear nature of (5). 

It was realized soon that what is needed is an approximation to the probability 
density function (pdf) of £, and that the mean value of any quantity which can be 
expressed as a function of £ is (Lin & O’Brien 1974, Bilger 1976) 

(/) = J KOpiOdi, (6) 

where p(£) is the pdf. Numerous experimental (LaRue k Libby 1974, Anselmet k 
Antonia 1978, Breidenthal 1981, Mungal k Dimotakis 1984, Koochesfahani k Di- 
motakis 1986), theoretical or numerical (Eswaran k Pope 1988, Pumir 1994, Holzer 
k Siggia 1994), and modeling (Kollman k Janicka 1982, Broadwell k Breiden- 
thal 1982) efforts have been undertaken to understand the properties of the pdf of 
passively mixed scalars. 

Of particular interest in this report is the /3-pdf model of (Cook k Riley 1994), 
in which the form of the scalar pdf is modeled as a function of its mean value and 
standard deviation and, especially, its use as a sub-grid model for large eddy sim- 
ulations (LES). Large eddy simulation has proven to be a powerful technique for 
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the computation of complex flows and good results have been obtained in the com- 
putation of filtered scalar mean values (Lesieur &: Rogallo 1989, Moin et al. 1991). 
We will show below that the subgrid fluctuation intensity can also be estimated 
with good accuracy. In this report we will assume that exact filtered mean scalar 
values can be computed by some LES scheme, but we will obtain them by filtering 
direct numerical simulation fields. The /?-pdf model has been tested in this way 
for isotropic turbulent flow at relatively low Reynolds numbers in (Cook Sc Riley 
1994, Reveillon Sc Vervisch 1996). We will test it here in the more realistic case of 
a mixing layer at medium Reynolds numbers (Rogers Sc Moser 1994). 

At issue is the question of large-scale turbulent inter mi ttency, which is the pres- 
ence of essentially laminar pockets in an otherwise turbulent flow, and whether the 
same subgrid mixing model can be used in homogeneous turbulence and in the pre- 
sumed interface between turbulent and laminar regions. A lot of effort has gone 
into modeling such intermittency effects (Libby 1975, Dopazo 1977, Ivollman 1984, 
Pope 1985, Pope Sc Correa 1988) but, if it is really a large scale effect, LES should 
be able to resolve it without resorting to modeling. The main difference between 
homogeneous flows and the mixing layer is that, while large-scale intermittency is 
rare in the former, it is prevalent in the latter. 

The simulation experiments are described in the next section. The results of 
applying the /?-pdf model to the prediction of different quantities in infinitely fast 
chemistry are presented in §3. We then extend the model to finite rate chemistry 
in the flamelet limit, and introduce some results on the joint pdf of the scalar and 
the scalar gradients, followed by discussion and conclusions. 

2. Numerical experiments 

The two flow fields used in this report are taken from the simulations in (Rogers 
Sc Moser 1994), where they are described in detail. Briefly, they are direct simu- 
lations of three-dimensional, temporally growing mixing layers, spatially periodic 
in the streamwise and spanwise directions, with initial conditions which represent 
turbulent boundary layers. The flow fields chosen are those in Figs. 18. a and 18. c of 
that paper, at which time the momentum thickness, 6, of the layers has grown by 
factors of 2.47 and 2.94 respectively from the initial conditions, and the streamwise 
integral scale has increased by a factor of about four. The energy spectra have a 
short power-law range with an exponent close to —5/3, and the layers are grow- 
ing self-similarly. Both layers appear to be slightly beyond the “mixing transition” 
identified in (Konrad 1976, Breidenthal 1981). The Reynolds numbers based on 
the instantaneous momentum thickness are 1980 and 2350, and correspond to lon- 
gitudinal microscale Reynolds numbers Re \ = 127 and 214 at the central plane of 
the layer. Note that both flows are quite intermittent, especially the second one, 
and that these values would probably change if they were conditioned only to the 
turbulent fluid. The ratio between vorticity and momentum thickness is about 4.85 
in both cases. 

The computational boxes are, in both cases, 125 x 31.25 initial momentum thick- 
ness and contain five or six large spanwise structures at the times chosen for our 
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FIGURE 1. Scalar intermittency across the two mixing layers used in the text, 

defined as the fraction of fluid for which £ £ (0.02, 0.98). : unforced case. 

: forced. 

experiments. A passive scalar is introduced at the initial condition with a laminar 
profile of thickness similar to that of the vorticity layer and a range £ £ (0, 1). 
Although the original simulations were spectral, using a mesh that was neither 
uniform nor isotropic, they were spectrally interpolated to physical variables on a 
uniform isotropic mesh for the purpose of our experiments. This implies some re- 
duction in the resolution, which is dictated by the least resolved direction. Thus 
while the original computations were carried using 512 x 180 x 128 and 384 x 96 x 96 
spectral modes, the interpolated fields contain 512 x 128 x 128 and 384 x 120 x 128 
points. The pitch of the interpolated grids is Ax/0 ^0.1 in both cases, although 
the original grids are finer by about a factor of two, especially at the central plane 
and in the transverse, y, direction. All lengths in this report are normalized with 
the instantaneous momentum thickness of the layers. In terms of the Kolmogorov 
scale at the center of the layer, 0/rj — 67 and 72 respectively, and the resolution of 
our interpolated grid is about 7 r) in both cases. 

It was found in (Rogers & Moser 1994) that the structure of the layer depends dur- 
ing the whole simulation on the initial conditions, corresponding to similar long term 
effects during the initial development of experimental layers. Different amounts of 
initial perturbations were introduced in the simulations to mimic this effect. Our 
two flow fields correspond to two extreme cases in the amount of two-dimensional 
perturbations applied to the initial conditions. In the first one, which will be re- 
ferred from now on as the “unforced” case, the initial conditions were synthesized 
from two turbulent boundary layers without modification. At the time of our ex- 
periment, both the vorticity and the scalar field are fairly disorganised with weak 
spanwise coherent structures, and there is very little fresh fluid at the center of 
the layer. In contrast, the second case was initialized by amplifying the spanwise- 
coherent modes of the initial boundary layers by a factor of 20, and the resulting 
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two-dimensional forcing gives rise to clear spanwise rollers with fresh fluid from one 
or the other stream present across the layer. This is clear in Fig. 1, which presents 
the fraction of mixed fluid in both cases, arbitrarily defined as f € (0.02, 0.98). 
Not only is the mixed fraction higher in the unforced case, but the presence of a 
few larger structures in the forced one results in insufficient statistics which are not 
symmetric with respect to the central plane. The statistics for the unforced case 
are symmetric. 

Of the two cases, the forced one is the hardest to compute because of the larger 
intermittency. Most of the results given below are for this case. The corresponding 
ones for the unforced case are at least as good, and usually better. 

We will generally compare mean quantities, denoted by {•), which are averaged 
over whole x — z planes. Occasionally the averages will be extended to slabs of 
the mixing layer, in which case the limits of the transverse coordinate y will be 
given explicitly. In our simulations of LES we define our basic filtering operation 
as a box filter in physical space. Quantities are averaged over a cubical box of 
contiguous grid points of side h = nAx, and assigned to the center of the box. 
This operation will be denoted by an overbar. Other filtering kernels have been 
used by other investigators, and it is not clear which is the best choice to mimic 
the projection operation implicit in a discrete grid, but our choice seems natural 
for finite differences or finite volumes codes, and has the advantage of providing a 
simple definition for subgrid statistics. 

Equation (6) extends trivially to filtered quantities, but the pdf has then to be 
taken to refer only to the interior of the filter box. Thus for a filter of width h we 
can define a subgrid mean 

£(x) = /T 3 J £(x-x')d 3 x' = Jtph((i;x) d£, (7) 

and a variance as 

= ( 8 ) 

All quantities are functions of y and, in addition, filtered quantities are also functions 
of the homogeneous coordinates x and 2 . To increase the number of data points 
available for the statistics, filtered quantities are computed at all grid points, even 
if they are only strictly independent over a coarser grid of pitch h. Plane averages 
are then computed for these filtered quantities and used to generate filtered profiles, 
which satisfy 

fh / 2 

(f) = h~' {/) dy = (/). (9) 

J-h/2 

Note that we can combine (6) and (9) to generate a “filtered” pdf for 

fhj 2 

y) = fc -1 / p{£,y - y')dy', 

Jh / 2 


( 10 ) 
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FIGURE 2. Mixture fraction pdfs from LES without subgrid modeling. Forced 

layer, (a) y/6 = 0. (b) y/9 - 3. : no filter. : h/9 = 0.44. : 0.66. 

: 0.88. Vertical arrows mark the mean value of £ for each plane. 

such that (6) generalizes to 


(/) = J mm<x- (ii) 

The error of any approximation to the mean profiles depends on our success in 
approximating p from our local models for p h . In most of our experiments the filter 
width will be small enough with respect to the width of the layer that we will be 
able to neglect the difference between laterally filtered and unfiltered pdfs. 

3. Fast chemistry 

3.1 No subgrid model 

It should be clear from the discussion in the last section that the aim of any 
approximation should be to reproduce p(£) as closely as possible. In RANS compu- 
tations, all the available information is the mean value of the mixture fraction over 
a plane and perhaps some of its statistical moments. Unless some model is applied 
for the form of the pdf, the implied representation is a delta function p = <$(£ — {£}), 
and is known to be poor. 

In LES we have some hope of avoiding subgrid modeling, since the grid elements 
are small parts of the flow in which the fluid may be assumed to be mixed and 
well represented by its mean. Large intermittent unmixed regions are hopefully 
contained in individual grid elements. In this approximation 

p*(0 = «(*-?), 7 (0 = /(€)- (12) 

In practice the filtered grid values are treated like real points and used to compile 
statistics. 
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FIGURE 3. (a) Sketch of the two filters for the estimation of the subgrid vari- 

ance. (b) Joint pdf of the band-passed mixture fraction fluctuation L and its true 
subgrid value £' h . The dashed line is (15) with spectral slope — 5/3. Isolines are 
logarithmically spaced by half an order of magnitude. Forced layer, y/6 = (—2, 2). 

The approximation (12) is tested directly in Fig. 2 for the pdfs of the mixture 
fraction in two planes of the forced mixing layer. Each figure contains the pdfs 
resulting from several different filters, compared to the real one. The widths of 
the filters axe of the order of the momentum thickness (30-50 Kolmogorov lengths), 
and correspond to grids of 10-20 points across the layer (Fig. 1). Even with these 
relatively coarse grids it is interesting that the approximation of the pdf is already a 
large improvement over the delta function of the global mean, and that the general 
shape of the pdf is recovered. Product mass fraction profiles obtained from using 
these pdfs in (11) have errors of the order of 20%. 

Nevertheless there are clear differences between the true and the approximate 
pdfs. A sizable percentage of the pure fluid that should be associated to the delta 
functions at £ = 0 and £ = 1 has been aliased as mixed fluid into the central peak. 
This is especially evident deep into the layer (Fig. 2. a) where the unmixed regions 
are presumably of small size and are almost completely obliterated by the filter. 

3.2 The Beta subgrid model 

To improve the approximation in the previous section it was noted in (Cook & 
Riley 1994) that, if the subgrid variance (8) were known at each grid point, it 
should be possible to make a reasonable guess as to the form of the subgrid pdf, 
Ph(£] £, £JJ, and to obtain a better estimate of the true pdf in terms of the joint 
pdf of those two subgrid variables 

m = J p(i &)p*(e; i (13) 

In the particular model proposed in that paper, the subgrid pdf is represented as a 
Beta distribution, ph ~ £ a_1 (l — £) b ~ x , and the two exponents a and b are computed 
at each point from the values of £ and £^ . 
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FIGURE 4. (a) Proportionality constant between band-passed and subgrid fluc- 

tuations, versus spectral slope, (b) Compensated scalar longitudinal spectra at the 

central planes (arbitrary units). : unforced layer. : forced. Vertical 

dotted lines are kh = 7r for the two filters used in Fig. 3.b. 

This correction needs the subgrid scalar variance, fj,, which is generally not given 
by the LES equations, but the same paper suggests that it may be obtained by 
a similarity argument from the behavior of the scalar at scales close to grid filter. 
Consider in our implementation two filter levels (Fig. 3. a). The first one is the grid 
filter of width h, which is represented by the dashed squares. The test filter £ is 
formed by averaging a 2 3 cubic box of contiguous, non-overlapping, grid values. In 
this implementation 

t = l (14) 

and we can define the subgrid variance at the test level 

tfh=T 2 ~f- (15) 

Neither (15) nor (8) are known, but they can be combined to give a band-passed 
“Leonard” term which, using (14), can be written as 

L 2 = tf h -tf=? -t- (16) 

The right-hand side involves only filtered quantities, and can be computed as the 
standard deviation of £ within the box defining the test filter. The similarity as- 
sumption is that 

^ = clL, (17) 

and is seen to be reasonable in Fig. 3.b, where it is tested for the central part of 
the forced layer. 

The proportionality constant can be estimated by assuming a form for the scalar 
spectrum, S(k ) ~ k~^ . The subgrid variance is obtained by filtering the spectrum 
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through the transfer function of the filter, which has the form F(kh). The result is 
that £'h ~ h B ~ l and, from (17) 


ci — (2 /?_1 — l) -1 / 2 . (18) 

This quantity, which assumes an infinite Reynolds number in that it integrates the 
spectrum to k -+ oo, is plotted as a function of the spectral slope in Fig. 4. a. For the 
Kolmogorov slope (3 = 5/3 it has a value ci = 1.305, which is the one used for the 
dashed line in Fig. 3.b, and represents the data well. In reality it is known that scalar 
spectra have slopes which are somewhat lower than 5/3 for Reynolds numbers in 
the range of our experiments (Sreenivasan 1996). This would imply proportionality 
constants somewhat higher than our value, but this effect is partly compensated by 
the presence of a Kolmogorov cutoff in the spectrum, which would lead to a lower 
value of ci. Figure 4.b shows our scalar spectra and the position of our filters with 
respect to them. The fact that both effects compensate at our Reynolds number, 
and that they should vanish as Re — ► oo, suggests that the asymptotic value of ci 
is a reasonable approximation for most of the Reynolds numbers of interest in LES. 

In the two previous tests of the /3-pdf model, the proportionality constant ci 
was fitted to the data and found to be smaller than ours. Reveillon and Vervisch 
(1996) found ci = 0.5 for a filter ratio of two, while Cook and Riley (1994) found 
cl ~ 1 for h/h = 1.8, which would correspond to /3 ^ 2.15 according to (18), and to 
cl ~ 0.9 for h/h — 2. Both simulations, however, were carried at Reynolds numbers 
substantially lower than ours. Reveillon and Vervisch worked at Re a = 17, for which 
there is no inertial range and no self-similar spectrum, and where turbulence is still 
barely developed. Cook and Riley do not give their Reynolds number, but their 
filters are only 6 times larger than the Batchelor scale, which would be near the 
right-most points in the spectra in Fig. 4, and within the dissipative range. Neither 
experiment can therefore be expected to agree with an inertial range prediction. It 
is probably a general rule that, if LES models are to behave independently of the 
type of flow, they should only be used in well- developed turbulence with filters in 
the inertial range. 

The results of applying the (3 correction to the pdf in the previous section are 
shown in Fig. 5, where it is seen that error has decreased considerably with respect 
to Fig. 2 and, especially, that it is now relatively insensitive to the filter width. Note 
that the good behavior of the model is not only at the level of integral quantities, 
but at the detailed level of the pdf, implying that it should give good results for 
the average of any function of £ and not only for the mass fraction (5). This 
includes the approximation of the pdf at a particular value of f, which is useful, 
for example, in evaluating source terms located at the flame. The figure includes 
the Beta distributions corresponding to the global averages and (true) standard 
deviations at each plane, as would be used in RANS. 

Although they are not included in the figure, there is no appreciable difference 
between the LES results obtained using the true value of and those obtained 
using the estimation (17). 
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FIGURE 5. Mixture fraction pdfs from LES using the similarity /? subgrid model. 

Forced layer, (a) y/6 = 0. (b) y/9 = 3. : no filter. : h/6 = 0.55. 

: 0.77. : 0.99. -o- : Beta distributions using the global mean and 

standard deviation for each plane. 


It is curious that, when the RANS pdfs are used to compute the mean value of 
the mass fraction (5), the result is within a few percent of the actual one, but it 
is clear from Fig. 5 that this is due to compensating errors and that it cannot be 
extrapolated to other quantities. 


3.3 Mean profiles 

The results of using the approximate pdfs of the previous section to compute 
mean profiles of various quantities are presented in Fig. 6, in which the degree of 
difficulty increases from top to bottom and from the left to the right. Plots on the 
left of the page are computed for a stoichiometric mass fraction £ s = 1/2, for which 
the flame is roughly in the middle of the mixing layer. There the fluid is relatively 
well mixed, and the results should be comparable to those obtained in homogeneous 
turbulence. Those on the left of the page are for £ s — 1/9, which corresponds to 
global models of the H 2 -O 2 reaction. For this stoichiometry the flame is near the 
edge of the mixing layer, in the interface between mixed and unmixed fluid, and 
LES may be expected to have more problems. The first two plots are mass fraction 
profiles obtained from the relatively smooth function (5). Those in the middle are 
profiles of Yp 4 , which is proportional to the fourth power of the temperature, and 
would therefore be a rough model for radiative heat in a flow with a real, hot, flame. 
This function (see Fig. 7. a) is much sharper than Y p and is therefore sensitive to 
the local values of the pdf, in spite of which the errors in the mean profile are still 
small. The last two plots are the values of the pdf at a given £ s and are the most 
sensitive test of the three. They are also the ones for which the errors are larger, 
but it is remarkable that the general form of the profile is still captured and that 
the errors stay, at worse, of the order of 25%. 
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FIGURE 6. Average profiles for different functions in the forced layer, (a)-(b): 
Product mass fraction Y p . (c)-(d): “Radiation” source Y*. (e)-(f): Pdf(f s ). (a), 
(c) and (e) are for = 1/2. The other three are for, = 1/9. 
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FIGURE 7. (a) The three functions used for the profiles in Fig. 6. £* — 1/9. (b) 

Relative variation of the profile thickness with filter width. Forced layer. Lines with 

symbols are for = 1/9. Other lines are for = 1/2. : Y p . : Y* 

: Delta function with = 0.02. 

For any profile which vanishes at y = ±oo we can define a ‘‘thickness” 

/ OO 

f(y)dy, (19) 

-OO 

which is proportional to the total amount of the particular quantity contained in 
the layer, and which can be used to quantify the global error of the approximation. 
Note that this thickness is unchanged by the filtering, 6j — Sf. The results for the 
different profiles of Fig. 6 are presented in Fig. 7.b, where they have been normalized 
with their DNS values. The errors for h/8 < 1 vary from better than 5% for the 
product mass fraction, to about 15% for the pdf. They are, as expected, generally 
larger for flames near the interface than for those at the center of the mixed region. 

4. Finite rate effects 

If the speed of the chemical reaction is large but not infinite, it is still possible to 
treat the combustion problem as a perturbation of the Burke-Schuman limit that 
we have used up to now. In this “flamelet” regime the deviations from infinitely fast 
chemistry are confined to a thin region around the location of the stoichiometric 
mixture fraction, whose width is a function of the Damkohler number. 

Although the reaction zone is typically thin, there are cases in which the non- 
equilibrium effects are globally important. One such example is the H 2 -0 2 reaction, 
in which an intermediate species is the H radical which, even in small amounts, 
controls the global exothermic properties. A simplified scheme 

3 H 2 + 0 2 — ► 2 H 2 0 4- 2 H, (20. a) 


2 H -j- M — » H 2 -j- M, 


( 20 . 6 ) 
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FIGURE 8. Pdfs of the mixture fraction gradients conditioned on the mixture 

fraction and compiled at different locations across the forced layer. : £ = 0.2. 

; o.5. : 0.8. : 0.2. -o- : 0.5. The first three curves are for 

y/6 G (-1, 1), (0 = 0.46. The last two are for y/6 € (1, 3), (£) = 0.72. 

was analyzed in (Sanchez et al. 1995). The first reaction is very fast and can be 
described by a conserved mixture fraction £ and by an infinitely thin flame located 
at its stoichiometric isolevel £ s , while the second one is slower and leads to finite re- 
action rate corrections. The Damkohler number is defined a s D — (Kx 2 t r )~ 1 , where 
k is the diffusion coefficient of the scalars, \ = is the mixture fraction gradient 
at the stoichiometric level, and t r is a chemical time which is only a weak function 
of temperature. The combination k\ 2 is usually called the scalar dissipation. In 
our approximation the only variable is \ , which therefore controls the structure of 
the flame. Other reactions involving radicals are technologically important. For 
example, the NO x production in air is controlled by the temperature and by the 
concentration of the O radical. 

It turns out that both the thickness and the maximum concentration in the radical 
containing region are proportional to D , so that the total mass of H radical per 
unit flame area is proportional to 

m ~ D~ 2/3 ~ k 2/3 x 4/3 - (21) 

There is a chemical energy associated to this mass which leads to a lowering of 
the flame temperature. The dependence on a power of the gradients is common 
in many other examples of slightly out-of-equilibrium reactions (Williams 1985), 
although the exponents change from one case to another. Assume in general that 
m — K a x 2a - We can estimate the average mass fraction of radical by a procedure 
similar to that used to derive (6). If S is the area per unit volume we write 
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which we wish to transform into a probability integral. Introduce the joint pdf of \ 
and £ and define dn as the element of length normal to the flame, located at £ = £ s . 
We can define the volume element both in terms of the geometry and of the pdf 


dV = dndS = p 2 (x , £) d* d£, (23) 

from where, using that \ = d£/dn, it follows that dS — \p(\, £)dx, anc l 

(Y H ) = [ X^(.\)P2(X, 0<*X =p(6) / xm(x)p(x|6)dx- (24) 

*/£=£, Jo 

The new pdf which appears in the second part of this equation is the pdf of gradients 
conditioned to £ = £ a . Note that (24) has the same form as (6) for a function 


1*00 

f(0 = C x 6(^ - £,), C x = / xm(x)p(x\£s) dx, (25) 

JO 

which is a delta function at the location of the flame, with a prefactor which depends 
on a moment of the conditional gradient pdf. 

Formulas of this type have been known for a long time in the context of nonequi- 
librium chemistry (Bilger 1976), and the joint pdf of the gradients and of a mixing 
scalar has been the subject of intensive study. There is general consensus that the 
unconditional gradient pdf is approximately log-normal (Kerstein & Ashurst 1984, 
Anselmet & Antonia 1985, Eswaran & Pope 1987, Pumir 1994, Holzer & Siggia 
1994), a form for which there is incomplete theoretical support (Gurvich & Yaglom 
1967, Meyers Sz O’Brien 1981) but which seems to be only an approximation to 
the real one. There is less consensus on the conditional pdf and, in particular, on 
whether the conditional variance of the gradients is correlated to the value of £ a or 
to the local turbulent dissipation e. 

We have checked conditional gradient pdfs for our two shear layers and the results 
are shown in Figs. 8 to 10. It is seen in the first figure that the form of the pdf 
is fairly independent of both £ s and of the location in the flow, when each pdf is 
referred to its own standard deviation x # * If a ls° has a general log-normal shape, 
but is not really log-normal. Note that the figure includes pdfs conditioned on 
values of £ s close to zero, but compiled at locations at which the mean value of £ is 
close to one. 

Figure 9 presents a two-dimensional map of the conditional x 7 (0> as a function 
of £ and of the location across the layer. It was found that the map was more 
uniform if the conditional was normalized with the unconditional standard 
deviation at the center of the layer, than when the normalization was done with 
the unconditional a f the particular y location. The first choice is used in the 
figure. It is seen that the dissipation has a central plateau, in which xVXo ~ 
but becomes larger near the edges of the layer, and vanishes at £ = 0 and £ = 1. 
The latter is an obvious property of laminar unmixed fluid and will be discussed 
below. The decline is, on the other hand, quite local and only happens when £ is 
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FIGURE 9. (a) Standard deviation of the scalar gradient magnitude as a function 

of y/6, conditioned on £ and normalized with the unconditional \q at the central 
plane. Isolines are \'/\o = 0.6 (0.2) 1.2, and shaded area is \7\o > 1-4. Forced 
layer, (b) Pdf of £ for the same flow. Each horizontal section represents the pdf 
over one plane. Isolines p(£) = 0.1 (0.2) 2.1. Isolines alternate line style for clarity, 
and shaded regions correspond in both figures. 

within 10% of the unmixed fluid. It should therefore not be important unless the 
stoichiometric is very close to 0 or 1. The rise near the edges of the layer is 
real, but it corresponds to combinations of mixture fraction and location which are 
relatively improbable, as can be seen in Fig. 9.b, in which the areas of high scalar 
dissipation have been overlaid on a two dimensional map of the mixture fraction 
pdf. It is clear that they correspond to events whose probability is mostly below 
10% and which will not have a large weight in (24). 

In Fig. 10 we have presented cross-stream profiles of \'(0/\o f° r various values 
of £. These are essentially vertical cross-sections of Fig. 9. a, but they have been 
included to give some quantitative information on the magnitude of the deviations 
of the scalar dissipation from its unconditional value, and to present data from the 
unforced layer. As in the previous figure it is seen that the scalar dissipation in the 
central part of the layer, where the fluid is well mixed, is more or less constant and 
equal to its unconditioned maximum value, but that the gradients conditioned on 
mixtures fractions close to the free stream values and all the gradients near the edges 
of the layer have standard deviations that may differ from the global maximum by 
almost a factor of two. They also have a characteristic parabolic shape. Most 
of these high deviations occur at places at which the absolute probability of £ is 
small, as seen on the figures on the right hand side of 10, which are conditioned on 
p(£) > 0.1, and they will only have a small effect on (24), but the effect is real and 
begs, at least, for some theoretical explanation. 

It is also interesting to note that the general magnitude of the gradients is low. 
The value of \q0 is 1.5 for the unforced case and 2.0 for the forced one, so that 
the peak of the distributions in Fig. 8 is for gradients of the order of those of the 
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FIGURE 10. Variation of the conditional standard deviation of the gradients across 

the layer, as a function of £. : £=0.1. : 0.3, : 0.5. : 0.7. 

-o-: 0.9. (a) and (b) are for the unforced layer; (c) and (d) for the forced one. 
(a) and (c) present full profiles, but (b) and (d) are only for those points in which 
p(0> 0.1. 

mean £ profile. Since high gradients lower the Damkohler number and may lead to 
extinction of the flame, the probability of local extinction for a given flow can be 
read from these distributions. 

In LES computations it might be harder to estimate the value of \q than that 
of the subgrid scalar fluctuation. The problem is that while the spectrum of £ 
decreases with wavenumber approximately like S(k) ~ that for the gradient 

increases as k 2 S(k) ~ k ly/3 . Thus, while most of the contribution to £' comes 
from the resolved large scales, comes mostly from the unresolved small ones. In 
terms of the standard LES or modeling equations, £' 2 is equivalent to the subgrid 
energy, while \' 2 is equivalent to the subgrid dissipation. Conservation equations 
and closures for the subgrid dissipation have been written among others by Newman, 
Launder Lumley (1981) and Elgobashi &; Launder (1983). 

The “engineering" consequences of the errors due to the gradient pdfs are sum- 
marized in Fig. 11. Assume that we are interested in computing the total amount of 
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FIGURE 11. Computed scalar dissipation thickness as defined in text, compared 
to its DNS value, h/0 % 0.9. : unforced layer. : forced. 


H radical in the example (20), and that we use as a subgrid model for the gradients 
an average pdf taken from Fig. 8, and a representative unconditional \o estimated 
for the center of the layer. Using (24) and integrating as in (19), we obtain a “rad- 
ical” thickness which is the integral across the layer of the a f 1 moment of the 
conditional pdf of the gradients, weighted with the pdf of the stoichiometric mixture 
fraction. This can be compared with the result of using the true distributions, and 
gives a global error due to the simplified assumption on the gradients. The result 
depends on the exponent a, but it is especially simple in the case of a = 1 since 
the a -f 1 moment is then proportional to the scalar dissipation, and the integral 
can be obtained directly from the data in Fig. 9. In this case the approximation 
is equivalent to taking everywhere Xo &s an approximation to x'(f 5 ). This normal- 
ized thickness is not very different from the results for a = 4/3, and is presented 
in Fig. 11. It is seen that, because the deviations from a universal distribution 
are mostly associated with places in which p(£ s ) is small, the final errors are still 
reasonable, especially for the unforced case, although they become 0(1) when the 
stoichimetric ratio approaches 0 or 1. 

The reason for this failure is clearly that we have not taken into account that 
gradients have to vanish when the scalar is very near its maximum or minimum 
value. Simple engineering models should be able to alleviate this problem, but they 
are beyond the purpose of this paper. It is, however, interesting to estimate the 
width of the region for which a correction needs to be applied, which either from 
Fig. 9 or 11 is in this case about « 0.1, but which can be related to the Reynolds 
number of the simulation. It follows from the form of the scalar spectrum that, for 
Sc = 0(1), most of x* is associated with scales of the order of the Kolmogorov 
length r), and that the scalar fluctuations at length t are Af* ~ where 
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FIGURE 12. Joint pdf for £ and the subgrid fluctuation Forced layer, h/6 = 
0.88. y/0 € (—2, 2). Isolines are p 2 (^ = 1(2)21. Dashed semicircle is the limit 

of possible (£, combinations. Lower dashed lines are the limits below which the 
Beta distribution looks like a single broadened spike. 

L< is an integral length. Since LJr) % R e (Tennekes & Lumley 1972), it follows 
that the scalar fluctuations which carry the gradients are of order 

A^^'Re- l/ \ (26) 

which in our case 0.03. As long as and 1 — £ 3 remain large with respect to 
this value, the small eddies should not be affected by the proximity to the level 
of the unmixed fluid, but if they are of the same order as (26), large gradients 
become impossible. This suggests that the width of the lateral bands in Fig. 11 
should decrease as the Reynolds number increases, but it would be interesting to 
get experimental confirmation of that estimate. 

5. Conclusions 

We have shown that relatively simple subgrid models for the pdf of a conserved 
scalar can be used to obtain useful engineering approximations to global quantities 
in LES simulations of reacting nonpremixed turbulent shear flows in the fast chem- 
istry limit. This is true even when the flow, in our case two different mixing layers, 
contains substantial intermittency. 

The magnitude of the approximation error varies from less than 5% for the total 
amount of generated products, to about 15% for the pdf of the scalar itself. When 
finite reaction rate corrections are introduced in the flamelet limit, the model has 
to be extended to the pdf of the scalar gradients, conditioned on the value of the 
scalar. We have shown that those pdfs have an approximately universal form and 
that they can be expressed in terms of a single parameter, the conditional scalar 
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dissipation, which varies little except at places in which the probability of finding 
mixed fluid is low. It is possible in those cases to obtain the global concentration of 
intermediate products (e.g. radicals) with errors which stay in the 20-30% range, 
except for reactions with stoichiometric mixture fractions very near those of the 
free streams. We have argued that the range of stoichiometric ratios for which the 
approximation fails should decrease with increasing Reynolds numbers. It should 
be clear, however, that even within this range the integrated quantities such as 
product concentration or radiation thickness are well predicted (Fig. 7). 

The particular approximation used in our experiments is the /? model of (Cook 
k Riley 1994), but it is clear from the lack of correspondence between actual and 
assumed pdfs that other models might work as well. This also follows from similar 
observations of Cook k Riley in their paper, and is in contrast with the situation in 
RANS, in which it is known that good subgrid models have to be used for the as- 
sumed pdf if any but the simplest quantities are to be computed accurately (Fig. 5). 
It is important to understand the reason for this difference, which is essentially con- 
tained in Fig. 2, where the scalar pdfs are reasonably well approximated even in 
the absence of a subgrid model. This means that most of the scalar fluctuations are 
associated to scales which are resolved by the LES, even for coarse grids like the 
ones used here. All that is left for the model is to correct situations in which the 
subgrid fluctuation is strong enough that the use of the average as a representation 
of the pdf is no longer appropriate. 

The situation would still be hopeless if those fluctuations were large enough to 
allow for a considerable latitude in the choice of subgrid pdfs, but this is fortunately 
not the case. Consider the (£, £j [ ) plane in Fig. 12. It can be shown that there can 
be no points above the dashed semicircle, and that pdfs that fall on the semicircle 
must be formed exclusively by unmixed fluid with £ = 0 and £ = 1. In the same 
way, pdfs on the horizontal axis are single delta functions of uniform fluid with 
£ = £. Pdfs near that axis are roughly spread deltas, and those near the semicircle, 
spread bimodals. The border between the two cases varies for different models, but 
it is always near the two intermediate dashed lines in the figure, which correspond 
to the /?-model. Below those lines, the pdf are bells, and almost any model should 
be equivalent. Pdfs within the two crescents correspond to spread deltas near one 
or the other free stream, and are also easy to model. Pdfs in the high-fluctuation 
central part of the diagram are harder, and are likely to depend on more than two 
parameters. 

We have overlaid on the diagram a typical joint pdf for £ and £^, for a relatively 
wide filter in the intermittent “hard” flow but, even in this case, most of the mass 
of the distribution is associated with pdfs within the easily modeled part of the dia- 
gram. The Beta distributions form a flexible set of pdfs which interpolate smoothly 
between the different cases, and they provide a simple numerical tool to evaluate 
the necessary integrals. This explains their practical success, but the reason why 
the approximation works lies in the small value of the subgrid standard deviations 
in Fig. 12. Since we have seen in Fig. 3 that these deviations can be estimated from 
large-scale quantities using an infinite Reynolds assumption on the upper limit of 
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the spectrum, it is unlikely that much higher values might be found in other flows. 

The small values of the fluctuations are also the reason why our relatively crude 
estimation of ^ works so well. Even large errors in this estimation have relatively 
small effects on the final results, and some experiments in which the estimated 
subgrid fluctuations where systematically increased or decreased by 20% did not 
show any appreciable differences with the results shown here. 

The convergence of the gradients also needs some discussion. It appears at first 
sight that, since the dominant contribution to \ n comes from the high end of 
the spectrum, the estimates for this quantity would depend of the value of the 
Batchelor scale, and would diverge at high Reynolds numbers. What is needed in 
(21), however, is not but the scalar dissipation k\ 72 , and it is easy to see that, 
for a k ~ 5 / 3 spectrum, this quantity is independent of Reynolds number. 
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